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We compute the quark number susceptibilities in two flavor QCD for staggered fermions by adding 
the chemical potential as a Lagrange multiplier for the point-split number density term. Since lesser 
number of quark propagators are required at any order, this method leads to faster computations. 
We propose a subtraction procedure to remove the inherent undesired lattice terms and check that 
it works well by comparing our results with the existing ones where the elimination of these terms is 
. analytically guaranteed. We also show that the ratios of susceptibilities are robust, opening a door 

_H ' f° r better estimates of location of the QCD critical point through the computation of the tenth and 

twelfth order baryon number susceptibilities without significant additional computational overload. 
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I. INTRODUCTION 



Quantum Chromo Dynamics (QCD), the theory of strong interactions, may have a critical point in the temperature 
(T)-baryon number density (or the baryonic chemical potential [Ib) depending on the number of light flavors it has. 
The search for this critical point is one of the major experimental goals of many heavy ion collider experiments world 
wide. In particular, the STAR experiment at the Brookhaven National Laboratory has begun an extensive Beam 
Energy Scan(BES) program [l|. It aims to scan the QCD phase diagram for baryon chemical potential between 20 to 
i— i, around 400 MeV, looking for the signatures of the presence of the critical point. FAIR at GSI, Darmstadt, Germany 
and NICA at Dubna, Russia are proposed to be operational in future with a key objective as the search for and 
the study of, a QCD critical point. It, therefore, appears natural to have a first principles theoretical exploration of 
the QCD critical point, with a possible prediction for its location which could serve as a useful reference for these 
large scale experiments. The physics near the critical point is essentially non-perturbative. Lattice gauge theory is 
. the most successful non-perturbative tool which can be used for such an exercise. Fluctuations of conserved charges 
' like the baryon number or strangeness are sensitive indicators of the existence of singularities in the phase diagram 
. like the critical point 0, H[. It is expected that the baryon number susceptibility would diverge at the critical point. 
t—( ' On a finite lattice, however, this quantity will show a peak at the critical point which would get sharper as the 
7— I \ lattice volume is increased. The direct computation of the susceptibility at finite baryon density is difficult due to 
~^ • the infamous "fermion sign problem" . One of the techniques to circumvent the sign problem at finite density is to 
compute the baryon number susceptibility as a Taylor series expansion in chemical potential near zero [4j-[6|. The 
radius of convergence of the series should yield an estimate of the location of the critical point 0j- For the precise 
estimation of the radius of convergence, one needs to compute ratios of as many higher orders of baryon number 
susceptibilities as possible. Current state of the art is the eighth order susceptibility on the lattice 8] . Higher order 
terms are important in determining the critical point, and extending to higher orders in Taylor series is therefore 
desirable although the explosion in the CPU time required is severely constraining. 

There are two major issues that have to be addressed for the efficient computations of higher order quark number 
susceptibilities(QNS). Firstly, using the standard method of introducing chemical potential on the lattice @-H2i the 
computation of QNS beyond eighth order with reasonably good precision becomes rather expensive. At each order, 
the fermion matrix inversions account for the maximum time involved in computing the QNS. As shown in Q, twenty 
inversions are necessary for computing the eighth order susceptibility; it would increase to about forty for the tenth 
order. For computing higher orders of QNS on the lattice, the number of matrix inversions would thus increase 
drastically thereby increasing the computation time. Secondly the current estimates of the susceptibilities beyond 
fourth order on the lattice tend to become statistically demanding as they are rather noisy. This is due to the fact 
that there are delicate cancellations in the expressions of QNS between different terms. Moreover, the number of 
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such terms itself increases progressively with the order of the susceptibility. Each term has to be computed with 
appropriate precision in order to ensure the cancellation of is free of computational artifacts. It would be desirable to 
reduce the number of such cancellations. 

In this paper we attempt to address both these issues. We use a staggered fermion matrix in which the chemical 
potential, fi, enters as a Lagrange multiplier, multiplying the point-split conserved number density term. It alleviates 
the problems mentioned above, as shown in [13]. For computing the eighth order QNS, only eight fermion matrix 
inversions would be needed [l3j as compared to the twenty required in the conventional method Q . In the early years 
of lattice computations this method of introducing fj, linearly was discarded because the second order QNS computed 
using this term leads to a 1/a 2 lattice term that diverges in the limit of vanishing lattice spacing, a. Modifications 
of the lattice action eliminated this term exactly for the free theory. In addition, they also ensured the correct 

Fermi surface on the lattice. These modifications lead to the difficulties in computing the higher order QNS though. It 
therefore appears a worthy attempt to focus back on the linear in fi action and devise ways to obtain physical answers. 
A method of successful elimination of the divergence maybe relevant in another context as well. The recently proposed 
overlap operator at finite density [bH ] with exact chiral symmetry on the lattice has such a linear dependence in /i. 
Indeed, the suggestions of [9KL2I] of divergence removal do not seem to work in that case. 

Our paper is organized as follows. In Section [TT] we briefly review the basic formulae of the various susceptibilities 
and show the differences in the two methods. In Section IIIII we suggest a procedure to remove the lattice artifacts 
in the QNS computed with the linear term in fj,. Free fermion results are shown to yield the correct continuum limit 
with it. We then compare our results for full QCD with those in the standard exponential method We show 
that our results are consistent with the existing results for different orders of QNS and the ratios of baryon number 
susceptibilities for a wide range of temperatures. A summary of our results along with possible sources of errors and 
their refinement is given in the last section. 



II. FORMALISM 

The quark number susceptibilities(QNS) for two flavor QCD are defined as, 

T d l+J \nZ(T,fj, u ,n d , ,m u ,m d ) 

Xij{»u,Vd = 77 t. (1) 

V dfj,i/j,j 

where Z is the QCD partition function, 

Z(T,n u ,Vd) = J VUe~ SG DetDiDetD}. (2) 

Sg is the action for the gluon fields. We use the standard plaquette action for Sq- D i} for each i = u,d, is the 
staggered fermion matrix in presence of finite quark chemical potentials. Several choices of the Dirac matrix are 
possible on the lattice at finite density. As in the continuum, the chemical potential can be introduced as a Lagrange 
multiplier corresponding to the conserved number density on the lattice in the point split form. This term is added 
to the standard staggered fermion matrix to yield the D(fi) which we use in this work: 
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= [^ U i( x ) S x, y -i - 1 l l U}{y)5 xy+l - (1 - ^a)r ti ul{y)5 x y+l + (1 + ^a)7 li Ui{x)S x y _ i + ma 8 XlV . (3) 

i=l 

r)i are the phase factors remnants of the gamma matrices as usual. Replacing (1 ± jia) by exp(±/^a), one obtains the 
popular method Q of introducing fj,. We shall compare our results with those obtained by employing the latter. 

In general, the chemical potentials [ii and fij, for quark flavors i and j respectively need not be the same. But 
since isospin is a good symmetry for QCD, we set the chemical potentials for up and down quarks to be the same, 
Mi* = Mrf = f 1 - Hence the baryon chemical potential is just /j>b = 3/i. The baryon number susceptibilities can then be 
expressed in terms of the quark number susceptibilities(QNS) Xij- For two flavor QCD, the expressions for baryon 
number susceptibility of n-th order, \b-> are 

X { b = 2 [* 40 + 2x31 + X22 l ' 

X { b = 47 [X60 + 4X51 + 7X42 + 4X33] , 
( 8) 1 

X B = gj [Xso + 6x7i + 16X62 + 26x53 + 15X44] ■ (4) 
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In this work we would be interested in computing the baryon number susceptibilities at [Xb = since these quantities 
appear in the Taylor series expansion of the second order baryon number susceptibility expressed in powers of fis > 



X2 0(p B ) _ X20(0) (4) /MS\ 2 (6) (Vb\ 4 (8) /MflV 

t 2 ~ r 2 +Xs UtJ +Xb \3t) +Xb \3t) 



(5) 



If a critical point exist then the baryon number susceptibility should diverge at that point in the continuum limit. 
The radius of convergence of this scries, should therefore determine the location of the critical point in the T-fig plane 
of the QCD phase diagram. It can be defined as 



r n = lim , 



(n+1) 

Xb 

->+3) ' 

Xb 



or r n = lim 



( 2 ) 
Xb 



1 l/n 



(n+2) 



(6) 



At the critical point, all the XiP are positive definite and the successive estimates of the radius of convergence agree 
with each other Q . It is clear from the above definitions of the radius of convergence that one needs to estimate more 
and more higher orders of the baryon number susceptibilities in order to locate it with precision and reliability. 

The QNS in Eq. |T]) can be written in terms of the trace of the derivatives of the Dirac matrix in Eq. ([3]) and 
its inverse at /i = 0. All the expressions for the QNS upto eighth order are given in the Appendix of Ref.[7[ as 
expectation values of O^ki..- These in turn are written in terms of the derivatives and inverse of D. We note that as 
a consequence of changing to our lattice Dirac operator linear in chemical potential, only the expressions for Oijki.. 
change, and indeed simplify a lot. This is because the second and higher order derivatives with respect to the chemical 
potential of the D in Eq. §S§ vanish. For example, the expression of X40 m °ur case in the notation of Ref. Q is still, 



X40 = | [<0iin + 60ii2 + 40 13 + 30 22 + 4 ) - 3(0n + 2 ) 2 ] 



with each such O n simply given by 



0„ = (-l)"- 1 ^ - 1)! TT(D^D'y 



In order to appreciate the difference, let us point out as an example 02 and 04 in our case are: 

i-ln'\2 „„J /n _ cti/t-i-1 nM 
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6Tr(D~ L D'y 



while in Ref. Q they are 
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04 
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TriD^D") and 

- ITTxliD^D'fD^D"] - 3~b:(D- l D"f - 4Tr(D~ 1 D'D^ 1 D"') 



Tr(D- l D' : 



(7) 



(8) 



(9) 



(10) 



From a comparison of the Eq. (|9]) and Eq. (p~0|) . one sees the absence of the second derivative term in the former. 
In fact, it arises due to those modifications which eliminate the free theory divergence. But then due to the same 
reasons Eq. (flU)) has four additional terms of alternating sign as compared to Eq. ©. This is generic. As the order 
n increases, the number of such terms in the expression of the O n increases. It should be noted that each trace in 
the expressions of QNS contains product of the inverse of the Dirac operator with the derivatives of Dirac operator 
with respect to /i. The matrix inversion is the most expensive computation on the lattice. If successive derivatives 
of Dirac operator are all finite, one has to compute more number of matrix inversions. Using the operator defined in 
Eq. only the first derivative of the Dirac operator is finite. Thus the price of additional terms with sign changes 
increases at the higher orders. Computationally, this implies more inverses of the matrix D and more precision with 
each term has to be computed in order to get O n with comparable accuracy for all n. It therefore seems a worthy 
effort to devise a scheme to remove the free theory divergence in another way, as we attempt in this work. In one 
can successfully remove the lattice artifacts that appear in the expressions of the lower order susceptibilities, then it 
can potentially open the door to compute the eighth and the higher derivatives with considerably less computational 
effort, enabling a better check on the radius of convergence estimate, as argued in [l3T |. 



III. RESULTS 



In this section we first begin by computing the susceptibilities of free fermions using the staggered operator in 
Eq. ([3]) in order to separate the artifacts in the second and fourth order QNS. In [l3[ we proposed to remove the 
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FIG. 1: The contour diagram for calculating the number density for free fermions at zero temperature on the lattice. P denotes 
the pole at isinh -1 ^7.,, ,, . 



undesired artifacts by estimating the zero temperature contribution on a symmetric lattice for each value of the gauge 
coupling, as is done for the pressure or energy density computation. This was motivated by the idea of keeping the 
cut-off effects to be the same as temporal lattice size increases. On the other hand, if one notes that the all the 
analytic methods of removal of the /i 2 -divergent terms were proven to be so only for the free theory, and that no 
new divergences appeared in the simulations of the interacting theory, one is lead to believe that a numerical method 
devised also for the free theory may work as well. Furthermore, since no new renormalizations are necessary at finite 
T and /is in perturbation theory, we expect the free theory artifacts to be the dominant ones. In this work we show 
that this method of subtraction gives results for all QNS which are in good agreement with the existing results for 
T > T c . For T < T c , this subtraction method appears to lead to some differences but these are small enough to be 
tolerated as differences in the finite parts of two schemes of removal of infinity. On finer lattices these will shrink, 
if this is indeed so. Moreover, even for T < T c , we find that the crucial ratios of susceptibilities and therefore, the 
radius of convergence estimates are not sensitive to it. 



A. Free theory 



The number density and the quark number susceptibilities(QNS) can be calculated analytically for the free fermions. 
We sketch how it is done, and discuss the results for number density for free fermions, obtained with the fermion 
matrix in Eq. ([3]), and the popular exponential form. We consider a lattice with N sites along each spatial direction 
and Nt sites along the temporal direction. The lattice spacing is taken to be a and therefore the volume is N 3 a 3 and 
the temperature of the system is T = 1/(Ntcl). From Eq. (JTJ, the expression for the number density on the lattice is 

(sin uj n + i/iacosujn) cosw„ \ 



N 3 N' 



t ' — ' / + (sino; n + cosa; n ) 2 

p,n 



N 3 N T ; 



^2 F(u n ,iia,p) 



(11) 



where / = (ma) 2 + sin 2 (api) + sin 2 (ap2) + sin 2 (a£>3). For the other form, /x appears only as (u) n — ilia) in place of 
uj n above. The expression can be evaluated by the usual trick of converting the sum over energy states to a contour 
integral. The zero temperature limit for a fixed lattice spacing corresponds to , Nt — >• oo. The energy eigenvalues 
then become continuous and lie in the range [— , — ]. The expression for the number density in this limit is, 



— Y 

N 3 ^ 



-iF(p, fia)& ( 6 — sinh 



Vf 



(l-Z^a 2 ) 1 /* 



duj sina;(cosa; + i/xasinw) 



_ w 2n /+(1-^ 2 



i sin 2 to 



(12) 



where F(p,fj,a) is the residue of the function F(uj — iff) in the positive half plane and tan0 = a/i. In terms of the 
contour diagram in the complex uj plane in Fig. ([T]), the original term is the line integral 3. Completion of the contour 
leads to the residue and since the contributions of the line integrals 2 and 4 cancel with each other, one is left with 
the contribution of the line integral 1. It gives rise to additional lattice artifacts in the expressions of number density. 
It vanishes for the exponential form due to u) — > —lu symmetry. Note, however, that the residue is present in both 
cases, and leads to terms higher order in a in each of them. More specifically, in the continuum limit of a — > 0, the 
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number density in the two cases are respectively 
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(14) 



Comparing the expressions in Eqs. (JT3]) and (fl4l . one notices that the free theory divergence which, as mentioned 
earlier, exists for the Dirac matrix in Eq. ([3]). It is eliminated, on the contrary, analytically in the latter case. There 
is an additional difference which matters in any physical comparison. The second term of Eq. (TT2l contributes to the 



/i 3 term in Eq. (|13[) as well, and reduces its value. On the other hand, such a Fermi surface violation is also cured in 
the usual method giving the correct coefficient for the /1 3 . The expressions are similar in form for higher order terms 
which affect the values of the sixth and higher order susceptibilities on a finite lattice. Indeed, the coefficient of the 
jj, 5 term is even larger in in Eq. (| 14|) . a trend which persists for even higher terms as well. Thus their approach in the 
continuum limit will be correspondingly slower. 

Since the QNS Xno(n = 0) of interest to us here are computed from the number density in Eq. (fT3")l . they too have 
these lattice artifacts in the form 0(a n ~ 4 ). As we saw above, for n > 6 they exist for both forms of the Dirac matrix, 
and are a bit less of a nuisance for the linear form than the exponential one. But, the X20 has an C(l/a 2 ) term which 
diverges in the continuum limit. It has to be subtracted before the continuum limit is taken. There is an additional 
O(a ) term in the expression for the fourth order susceptibility which gives an incorrect result. Clearly, removal of 
these artifacts is necessary. Since we can trace them to the second term in Eq. (fT2"j) . we propose to compute them 
directly numerically from it, and subtract. Thus taking one more derivative with respect to fxa, setting a/i = in the 
resultant expression, the subtraction term of X20 can be computed on a lattice of volume N 3 and infinite temporal 
extent by first analytically integrating over the energy eigenvalues along the temporal direction. Thus we numerically 
evaluate 



»<°>=-i^x;(i-\/^) 



(15) 



and 



X4o(0) = - 



3 V L 3 + 2/ 
4 N 3 ^ \ 1 



/ V 1 + / 



(16) 



as the vacuum subtractions at zero temperature. Eliminating these artifacts, the second and the fourth order suscep- 
tibilities are shown as a function of 1 /N T in Fig. ^ and compared with the corresponding results obtained using the 
exponential form. In these plots the aspect ratio N /Nt is fixed to be four as it yields already the thermodynamic 
limit. As Nt becomes larger, corresponding to the continuum limit, it is evident from the plots that our proposal for 
the artifact subtraction does ensure that these quantities indeed do have the correct continuum limit. The difference 
between the results for lattice sizes Nt < 10 are due to finite cut-off effects which are comparatively larger for the H-K 
method. This difference becomes more for the higher order QNS as shown in the plot of the sixth order susceptibility 
in Fig. ©. 



B. Interacting theory 

In this section we compute the nth order baryon number susceptibilities at vanishing \lb for two flavor QCD using 
staggered fermions using the operator given in Eq. © and the subtraction scheme explained above for the free 
theory. Note that we follow all the existing computations in this respect: the analytic cancellation of the divergence 
was shown only for the free theory and the same form used for the interacting theory. We used the same configurations 
used previously for estimating all susceptibilities up to the eighth order on a Nt = 6 lattice @. For details of the 
configurations and the scale setting, we refer the reader to the Ref. [8[ . The lattice size used was 24 3 x 6 and the pion 
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FIG. 2: The second (left panel) and the fourth order (right panel) susceptibilities for free fermions as a function l/iV|> for 
different methods of introducing \x. 
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FIG. 3: The sixth order susceptibility for free fermions as a function l/iV|> for different methods. 



mass was fixed at M w = 230McV, as there. The critical coupling was defined from the peak of the unrenormalizcd 
Polyakov loop susceptibility. As discussed in the Section UH the expressions of susceptibilities contain trace of fermion 
operator insertions. For computing each trace, 500 random vectors were used. This was also the optimum number 
of random vectors used in the earlier work Q for computation of the sixth and eighth order susceptibilities. Thus, 
essentially all the computational details were maintained to be the same as in |8|. Our expressions for the various 
O n are different and we use the subtraction scheme for %20 and %4o. Our results for QNS are compared with the 
corresponding results of Q, i.e., the exponential form and with analytic cancellation of the divergence for the free 
theory. Further, we also compute the ratios of our susceptibilities ans compare with the known results. 



C. Second order 



We compute the zero temperature value of X20 for free fermions using Eq. (|15[) . by numerically summing over the 
momentum modes on a 24 3 lattice since the spatial volume for our interacting case is 24 3 . We then subtract this 
quantity from the X20 values in the interacting theory computed using the O2 for our case. The left panel of Fig. 
(HI compares our results with those of [8|, labeled as 'GG'. At 0.92T C , the value of our baryon number susceptibility 
matches with the existing GG results within errors. This suggests that including our subtracted value the of (O2) 
is effectively the same. Since the difference in the two cases comes only from the Tr (M _1 Af") for Q and the 
subtraction term of Eq. (|15[) on a 24 3 lattice for our case, we can infer that their values are equal within the numerical 
precisions. Since one is evaluated in the interacting theory while the other in the free theory, this equality justifies 
our ansatz that interactions do not lead to further divergent terms, as expected also in perturbation theory. Once the 
divergent term is subtracted the difference between the results in the two methods are due to different cut-off effects. 
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FIG. 4: The second order flavor diagonal(left panel)x2o/T 2 and the corresponding one normalized by the Stefan-Boltzmann 
value on a Nt = 6 lattice (right panel) as a function of T/T c . 



The disagreement at the highest temperatures is consistent with the fact that there are different cut-off effects in the 
free theory limit i.e at asymptotically high temperatures. They would vanish only in the continuum limit or very 
large lattices. Again this can be easily verified by numerical computations in the free case. 

We therefore check the variation of the ratio X20/X20, SB computed as a function of temperature where the Stefan- 
Boltzmann(SB) value is computed on a Nt = 6 lattice for both the cases. We expect that if the cut-off effects in the 
interacting theory are not very different from the free theory for temperatures much larger than T c , this ratio should 
be independent of the subtraction procedure used. From the right panel of Fig. (Q| it is evident that the GG results 
of @ are consistent with ours for T > T c . At 1.92 T c , the deviation from the Stefan-Boltzmann value is about the 
same(~ 20%) for both the methods. 



D. Fourth order 



We follow the same subtraction procedure again for the fourth order susceptibility and compute the subtraction 
constant using the expression in Eq. (fT6|) on a 24 3 lattice. Comparing the X4B in the two cases in the left panel of 
Fig. ([5]), a good agreement is evident with the existing GG results. One does see, however, differences arising perhaps 
due to the finite size of cut-off a = 1/ 6T. In the right panel of Fig. ([5]) , the ratio of \ab and the corresponding Stefan 
Boltzmann value on a Nt = 6 lattice is shown as a function of T/T c . The deviations from the continuum Stefan 
Boltzmann value are consistent with the free theory results for T > 1.5 T c indicating that the ideal gas cut-off effects 
are likely to be more dominant than the interaction effects. At 1.92 T c , the ratio is away from unity by ~ 5% in our 
computation whereas the deviation is about ~ 15% for the GG results. At T = 1.21 T c , the deviations are larger 
than that for free massless fermions indicating larger effects due to interactions in the medium. Of course, continuum 
extrapolation are needed to make prediction on the nature of the QCD medium at T ~ 2T C , especially in view of the 
small differences observed. 



E. Sixth and higher order 



For sixth order and above, we follow the lesson learnt from the free theory and do not adapt any subtraction. The 
sixth and the eighth order baryon number susceptibilities are displayed in Fig. ([5]). We find again a good agreement 
with the GG results. At temperatures below T c the larger error bars in the GG results compared to our results 
arise due to the presence of terms with varying sign in the former. One has larger number of matrix inversions for 
computing such terms and one employs noisy estimators to determine them. From the plots it is also evident that the 
values of the higher order susceptibilities fall to zero rapidly for T > 1.5 T c . The regime between 1.2 — 1.4 T c may be 
still be sensitive to the critical fluctuations as seen from the large deviations of the sixth and eighth order fluctuations 
from the Stefan-Boltzmann values. In the continuum, the values of sixth and the higher order susceptibilities are zero 
for free fermions and are finite in QCD only due to the interactions. To understand how dominant are the interaction 
effects it is important to reduce lattice artifacts and perform a continuum extrapolation of such quantities. Using the 
Dirac operator in Eq. ([3]) these artifact effects are reduced as compared to the standard operator and also performing 
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FIG. 5: The fourth order baryon number susceptibility (left panel) and the same quantity normalized by the Stefan-Boltzmann 
value on Nt ~ 6 lattice (right panel) as a function of T/T c . 
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FIG. 6: The sixth(left panel) and eighth order(right panel) baryon number susceptibility as a function of T/T c . 



continuum extrapolation will be easier. 



IV. RATIOS OF SUSCEPTIBILITIES AND RADIUS OF CONVERGENCE 



The ratios of different QNS are sensitive indicators of the location of the critical point and the extent of the critical 
region. For massless QCD, the sixth and higher order susceptibilities should peak at T c with 0(4) critical exponents, 
implying indirectly the existence of the critical point. Since our input quark mass is quite large we should expect a 
crossover at T c . As mentioned in Section [TTl the ratios of susceptibilities are used to define the radius of convergence 
and hence are important for determining the location of the critical point. 

Since the ratios are independent of the volume of the system, it was suggested to use them in comparisons with 
data a heavy ion collision experiment where it is difficult to estimate the volume of the fireball formed. Also such 
observables can be used for comparing lattice results at finite volume with the experiments [lj|. The ratio of the 
fourth to the second order baryon number susceptibility is related to the product of kurtosis(«) and square of the 
variance(cr 2 ) of the baryon number. In the heavy ion collision experiments, the kg 2 of the net-proton number is 
measured as a function of the center of mass energy of the colliding heavy ion beams. A non-monotonic behavior of 
kg 2 as a function of beam energy is a possible signature of the existence of the critical point @ and is currently being 
probed at the RHIC Q. 

We compute the quantity K — na 2 and compare with the value of the same computed from the susceptibility 
values in the Ref. Q. Note that these are still at zero chemical potential, and are meant more for comparison of the 
two methods. The results are displayed in the left panel of Fig. (JT]). A very good agreement is observed between 
the results in these two different methods. This implies that the ratios of susceptibilities are not very sensitive to 
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FIG. 7: The kurtosis(left panel) and different radii of convergence(right panel) as a function of T/T c . 



the subtraction scheme used. We expect that the radius of convergence estimates too would not be very sensitive 
to the different subtraction schemes used. In order to verify this, different radii of convergence defined in Eq. © 
were computed at 0.92T C using our results and compared with the corresponding GG estimates. The result of the 
comparison is shown in the right panel of Fig. ([7]). As expected the radii of convergence computed using two different 
subtraction schemes are roughly in agreement with each other. This is promising as we can hope to now extend our 
analysis in the critical region with susceptibilities at tenth order and beyond this way. We hope to be able to check 
whether the current results on the location of the critical point are changed in a a significant way with such estimates 
stemming from QNS beyond the eighth order. 



V. CONCLUSIONS 



Ever since it was recognized that the lattice free theory leads to a~ 2 divergences in the continuum limit, mostly the 
exponential form for the chemical potential term was used, as it eliminates the divergence for the free theory. Such 
an exponential form, however, seems not possible for the overlap quarks which have the same chiral properties as the 
continuum. Demanding chiral invariance for them even at finite density, one obtains an overlap operator of the form 
linear in \i. Moreover, the linear form for any fermions, overlap or staggered, leads to simpler non-linear susceptibility 
expressions. Absence of canceling terms, and requirement of fewer Dirac matrix inversions make the linear form more 
suited for extension to higher order non-linear susceptibilities needed to estimate the radius of convergence, or the 
location of the QCD critical point . 

In this work we proposed a numerical scheme to remove the free theory artifacts inherent in the linear form. In 
particular, subtractions were introduced for the second and fourth order susceptibility which we suggest should be 
evaluated for the free theory at zero temperature but the same spatial volume. As one sees from our Figs. (|4I5[) , it 
works well. Indeed, the main difference from the popular exponential form stem from the corresponding free theory 
due to the coarse finite spacing used at temperatures T > 1.5T C . For the practically more relevant ratios, such as those 
in the left panel of Fig. ([7]) used as signature for critical point or in the right panel of Fig. ([7]), used for estimating 
the radius of convergence, our method works as well as the exponential form within errors. It would be interesting to 
extend this work to higher orders at the critical point, estimates of which are currently obtained with susceptibilities 
up to the eighth order. 
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